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1—1 Abstract 

I I In this paper the chaos persistence in a class of discontinuous dynamical systems of fractional- 

order is analyzed. To that end, the Initial Value Problem is first transformed, by using the Filippov 

7—i regularization [1], into a set-valued problem of fractional-order, then by Cellina's approximate 

^ selection theorem [21 [3], the problem is approximated into a single- valued fractional-order problem, 

which is numerically solved by using a numerical scheme proposed by Diethelm, Ford and Freed 

00 0]. Two typical examples of systems belonging to this class are analyzed and simulated. 

( Kkeyword Fractional derivative; discontinuous dynamical system; Filippov regularization; differ- 

ential inclusion; numerical method 

o 

,— I 1 Introduction 

. !__( Dynamical systems (d.s.), discontinuous with respect to the state variable, occur in a large number 

of problems, especially from mechanics (dry friction with stick and slip modes, impacts, oscillating 
systems with viscous damping, elasto-plasticity, alternatively forced vibrations, braking processes with 
locking phases), electrical engineering (electrical circuits and networks with switches, power electron- 
ics), automatic and optimal control theory, convex optimizations, game theory, non-smooth control 
systems synthesis, uncertain systems, walking machines, biological and physiological systems and ev- 
erywhere non-smooth characteristics are used to represent switches (see e.g. El [71 [U O US and some 
references therein). In other words, in the real word non-smoothness is a common phenomenon. The 
underlying mathematical models are generally described by a set of first-order differential equations 
with discontinuous components. Particularly, the discontinuity is due to the discontinuity of the state 
variable, of the associated vector field, of the Jacobian, or to higher-order discontinuity. 

In practical examples, the discontinuity appears because of piecewise continuous or switch-like 
functions (e.g. signum, absolute value, Heaviside function -also known as the "unit step function" - 
maximum function, etc.). 

On the other hand, there are many continuous d.s. which display fractional-order dynamics. These 
mathematical phenomena allow to describe a real object more accurately than do classical "integer" 
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dynamics. One of the main reasons to use the integer-order models was the absence of methods for 
fractional differential equations. 

The real objects or phenomena such as viscoelastic systems, dielectric polarization, electrochem- 
istry, percolation, material science, polymer modeling, theory of ultra-slow processes, electromagnetic 
waves, evolution of complex systems, financial processes, secure communication, etc. are generally 
fractional (see [ni[Tl[Il[Tl[ISl[I3[T71IIHl[Tn]). 

We consider in this paper a class of unified Filippov discontinuous d.s. of fractional-order and 
prove that the underlying Initial Value Problem (I.V.P.) admits solutions which can be numerically 
determined. We also prove numerically that in these systems of lower than third order chaos may 
appear (as it is known, in the case of integer order, according to the well known Poincar-Bendixon 
theorem, chaos appears only at systems of minimum order three). 

These systems are modeled by the following autonomous I.V.P. 

= g{x{t)) + A s{x{t)), x{0) ^xo,teI^ [0, oo), (1) 

with g : M" M", A ~ (^'i.j)„xn ^^^^ constant matrix, the vector function s : M" M" being 
composed of signum functions (the most encountered case), i.e. 



s{x) = 



/ sgn{xi) 



\ sgnia 



(2) 



q the fractional order, < q < 1 (although there are also some applications for 5 > 1, as is the case 
in many applications). 

The right-hand side depends on p G M, the control parameter. 

Hereafter, we impose the following assumption on the discontinuous right-hand functions 
(til):The right-hand side is piecewise continuous (continuous on a finite number of m open domains 
Di C i?", i — 1,2, ...,771, and has finite (possibly different) limits from different boundary points, 
i.e. bounded discontinuities), g being continuous (linear function in the great majority of practical 
examples) . The null set of the discontinuity (switch) points is M = R"\ U Di. 

At every Af point, the discontinuity surface is a hyperplane (plane for 77 = 3). 

The considered class of systems is autonomous. Thus, the I.V.P. ([I]) can be written as follows 

d^x 

— =g{x) + A s{x), x{0)=xo, tel. (3) 

To achieve our goal, we consider the I.V.P. ([s]) in a more familiar framework, such as differential 
inclusions and differential equations of fractional-order. 

The paper is organized as follows: In Section 2, the discontinuous I.V.P. ^ is transformed, via 
Filippov's regularization, into a set-valued I.V.P. and then into a continuous differential equation of 
fractional-order via Cellina's Theorem [H [3]. In Section 3 the numerical integration of differential 
equations of fractional-order, which will be used to integrate the obtained I.V.P. ([s]), is presented 
briefly. In Section 4 this algorithm is used to investigate the lowest-order q for which chaos vanishes 
in two representative chaotic systems. 



2 Continuous approximation of I.V.P. ([S]) 

One of the first questions of set-valued analysis is how to relate set-valued and single- valued functions, 
so as to avoid dealing with the set-valued functions. Because the considered examples in this paper 
are modeled by three-dimensional real I.V.P., the next results will be particularized in the Euclidean 
space E", some of them are valid even in more general (such as Hilbert) spaces. 
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First, let us consider the IVP ([s]) in the common case q = 1 (i.e. discontinuous LV.P. ) 

dx 

^g{x)+ A s{x), x{0)^Xq, tel. (4) 

at 

Due to the discontinuity of the right-hand side, the I.V.P. may not have any solution (see pj for the 
background on existence and uniqueness). 

To overcome this situation, Fihppov proposed the idea of restarting the I.V.P. as a set- valued I.V.P. 
via differential inclusion [T] 

dx 

— e F{x), x{0) — xq, for a. a. t £ I, (5) 

where F : E" =| M" is a convex set-valued vector function defined on the set of all subsets of M" (for 
the background on set-valued functions we refer e.g. to 3]). 

Let us consider the following discontinuous I.V.P. enjoying the assumptions HI 

dx 

— = fix), x{0) = xo, te I. 

One of the simplest definitions for F is the following 

F{x) = Pi Pi conv{f{x e M" :|| a; ||< e}\M), (6) 
e>o ^l(^I)=o 

F being the smallest closed and convex set containing all limit values of /, and /i the Lebesgue 
measure. 

On Di, where the function / is continuous, F(x) consists of one point, which coincides with the 
value of / at this point, F{x) — {f{x)}. Thus, on Di, the I.V.P. becomes a single-valued continuous 
problem dx/dt = f{x). At the discontinuity points, the set F{x) is a subset of M" given e.g. by (|6|. 

Remark 1 In the practical examples, e must he small enough, so that the motion of the physical 
system can he arhitrarily close to a certain solution of the differential inclusion. 

For example, the set-valued version of the usual sign function is defined by 

r {+1} a;>0, 
Sgn(x) = l [-1,1] x^Q, 
\ {-1} a;<0. 

By applying the Filippov regularization to the LV.P. Q, it becomes a set- valued I.V.P. 

dx 

— G g{x) + A S{x), x{Q) = Xq, for a.a. t £ I, (7) 



where S is the set-valued form of s i.e. S{x) — {Sgn{xi), Sgn{x2), . . . , Sgn{xn)Y 

Now, the obtained set-valued I.V.P. admits at least a generalized (Filippov) solution (for existence and 

uniqueness we refer to fl] and for numerical methods for differential inclusions [201 121j ) . 

Remark 2 It is easy to see that the set-valued function AS is upper semicontinuous with closed and 
convex values (see e.g. ISI P- 85 or p. 101). 

Next, following the Filippov procedure applied to the right-hand side of the I.V.P. (|3|, one obtains the 
following differential inclusion of fractional-order 

— e g{x) + A Six), x{0) = Xq, for a.a. t e I, (8) 
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Figure 1: Graph of the sigmoid function for different e values. 



Definition 1 f3f.l22f selection of a set-valued function F : M" =| M" is a single-valued function 
/ : M" — > M" satisfying 

h{x) e F(x), Vx e M". 

Convex set- valued functions admit approximate selections [3], i.e. single- valued functions : M" — > 
M" for which Graphih^) C Graph{F) + eB and h^{x) belongs to the convex hull of the image of F. 

Theorem 2 The set-valued function AS in the I. V.P. ^ admits a locally Lipschitzean approximate 
selection -.W — !■ M". 

Proof Because AS is upper semicontinuous with closed convex values (Remark [2|, according to Cel- 
lina's Theorem ( [5; p. 358, T. 9.2.1) it admits an approximate locally Lipschitz selection. 

The constructive proof of Cellina's theorem (see e.g. [2]) is a great advantage because it allows the 
construction of he ■ 

There are several strategies to choose selections for a general class of set-valued problem (see e.g. 

Following the preceding theorem, we can enounce the following main result 

Proposition 3 Let the I. V.P. with g continuous. Then, the L V.P. can he approximated by the 
following continuous LV.P. of fractional- order 

d'^x 

— = f{x), x{Q) ^xq, te I, 



with / : M" — > M" 



g{x) + he{x) for x € M , 
g{x) -f As{x) for x ^ M, 



where h^ is an approximate selection of the set-valued function AS. 
Proof 

The proof is obvious if we consider that first, the LV.P. ^ is transformed, via Filippov regular- 
ization, in the set-valued LV.P. ([s]). Next, Theorem [2] can be applied in the neighborhoods of the 
discontinuity points. 
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Remark 3 i) Obviously, s is a vector and in applications we have chosen the same value for all its 
components;; 

ii) As proved in 1231, this local approximation can be made even with smooth functions; 

Hi) The proof of the above theorem could also be done via maximal monotonicity of AS which ensures 

its closedness and convexity (see e.g. ]^ p. 14^1, Proposition 2). 

Because of the amount of necessary calculus for fractional numerical integration, we shall use a 
simple form for h^. 

To be precise, let us consider one component of As in the I.V.P. (|3|, the scalar function Sj{xj) = 
aijsgnixj). In a e— neighborhood oixj, we have chosen in our numerical experiments for the known 
sigmoid function, widely used in practical applications. Its scalar variant has the following form ( Fig. 

i 

To connect continuously in some e— neighborhood of Xj, we shall use the following form 

1 + e 

The influence of e is underlined in Fig. [5] d and Fig. [2] e, where the visible (for e = IE — 2) corners, 
typical for discontinuous d.s., become smooth if the width of the neighborhood enlarges. Therefore, 
for physical reasons e needs to be chosen relatively small. 

Conclusion The discontinuity hyperplanes of L V.P. can be replaced in their e -neighborhoods 
with continuous surfaces, via approximate selections. 

3 Numerical approach of discontinuous systems of fractional- 
order 

In this section we present summarily a numerical scheme for continuous differential equations of 
fractional-order which will be used for our I.V.P.([3]). 

3.1 Numerical integration of differential equations of fractional-order 

There are several definitions of the fractional differential operator. 

The Caputo fractional derivative of order g > 0, where q G (to — 1,to] for some to e N for 
a differential function x : [0, t] — > M, is the most common tool, widely used in real applications, 
allowing for interpretation in a physically meaningful way (see e.g. @], |15j and |24j ) 

1 } j;(")(s) 

D^x = — / — ds, 

T{m - q) J (t - s)«+i-™ 


where F is the known Gamma function. 

To be precise, let us consider in the sequel the continuous autonomous fractional I.V.P. 

D'ix = f{x), 0<t<T, x('=)(0) = 4"-\ ^^Q^ 

fc = 0, 1, . . . , TO — 1. 



Under the continuity condition of the function /, there is at least a solution (see e.g. [4] where the 
most general case of nonautonmous problems can be considered) . Because has an m-dimensional 
kernel (to being just the value rounded up to the nearest integer i.e. m — \q~\), m initial conditions 
need to be specified. Certainly, for the usual case < g < 1 we have to specify just one condition. 



5 



Table 1: Pseudo-code of DFF numerical scheme 



Xi[0] = a;io[0], i = 1,2, . . . , Ne {initial conditions} 

for j :— 1 to N do 

begin 

for i 1 to A^e do {predictor} 
hi i-1 

Pi ■- Xio[0] + — — — b[j - k] fr {xi [k\,..., XN, [k]) 

^ \1 + ^) k = 

for z := 1 to Ne do {corrector} 

X^ [j] ~ X, o[0] + ^^^.A(Pl, ■ • ■ ,PiVe) + ((i - 1)'+' - (j - 1 - q)f).f {^1 [0] , . . . , XN, [0]) 

i-1 

+ a[j -k] f {xi [k],..., XN, [k]) 

end 



The fractional differential operators are computationally time expensive as compared to their 
integer-order counterparts. Some numerical methods for ( |10[ ), essentially ad hoc methods for very 
specific types of differential equations are presented e.g. in ^1], pSl [251 [?7] . 

In agreement with the standard mathematical theory ([2H] Sect. 42), the initial conditions for the 



I.V.P. ( 10 1 should have the form 



However, because in practical applications, these values are frequently not available, and it may 
not even be clear what their physical meaning is, one can specify the initial conditions in the classical 
form as they are commonly used in initial value problems with integer-order equations [24] . 



The numerical scheme to integrate (10), which is used in this paper, is a predictor-corrector al- 
gorithm (called hereafter DFF), proposed by Diethlem, Ford and Freed in [3] together with error 
estimates and numerical examples. The scheme is a generalization of the classical multistep method 
Adams-Bashforth-Moulton and uses the Caputo derivative. 

Being of practical significance and helpful for solving a broad class of problems, the method has 
been constructed and analyzed for the fully general set of equations without any special assumptions, 
and is easy to implement on a computer 

Let us next assume that we are working on a uniform grid {i„ = nh : n — 0, 1, A^j with some 
integer TV with the step size h := T jN and the case of g G (0, 1]. 

The predictor formula for the value at the point i„_|_i,xP, is the fractional variant of the Adams- 
Bashforth method which for our particular case m = 1, has the expression 

1 " 

(i„+i) = + — ^ ^ ^j,n+i/ (a; (tj)) : 
while the corrector formula (the fractional variant of the one-step implicit Adams Moulton method) 

a;(t„+i) = a:o+ ^{q^l) ^^^'^ fa+i)) 



5^(^^^«i,n+i/(a^ (i,)). 



where F is the Gamma function, and a and h are the corrector and predictor weights respectively given 
by the following formula 
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- {n - q) {n + If if j = 0, 

-2(n-j + l)«+i ifl<J<n, 
1 ifj = n+l, 

and 

b,,n+i = j{{n + l-jr-{n-jf). (12) 

After a few modifications, designed to enhance the efficiency , DFF scheme has the form presented 
as pseudo-code in Table [l] for the case m = 1 and a set of autonomous differential equations of 
fractional-order q with the right-hand side 

where: 

Xio, i = 1,2, Ng are the initial conditions; 

x[-] array of A'^ -|- 1 real numbers containing the approximate solutions: Xk ~ x{kh), x being the 
exact solution; 

p the predicted values; 

a, 6, arrays of 3 x (N + 1) real numbers which, in this algorithm, are determined with the following 
variant 



for fc := 1 to iV do 

a[k] := ki -{k- If 
b[k] := [k + 1)"+^ - 2fc'=+i + (fc - 1)«+^ 
end 



The Gamma function T can be approximated via the so-called Lanczos approximation ^29) using 
the formula ( |http: / /www.rskey.org/gamma.htm[ ) 

6 

T{z) = + 5 5)^+o.5g-(.+5.5)^ (^3) 

n(^+*) 

for a complex variable z for which Re{z) > 0. The coefficients P are drawn in Table [2] 



Table 2: P coefficients for T function (13) 



^0 


= 75122.6331530 


Pi 


= 80916.6278952 


P2 


= 36308.2951477 


P3 


= 8687.24529705 


Pi 


= 1168.92649479 


P5 


= 83.8676043424 


Pe 


= 2.50662827511 



Remark 4 i) Another way to deal with fractional derivatives is the use of the frequency domain ap- 
proximation, based on the approximation of the fractional operators in the frequency domain. This 
technique, proposed by researchers on automatic control, is suitable for Simulink (Matlab) (see e.a. \3(K 
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a) As known, an attractor requires a relatively large integration time interval (for example for some of 
our experiments one needs N — 2EA^ 2.bEA), and an array of this size is hard to implement in any 
compilers. The best solution we found for our examples was to save x into a file. 



3.2 Numerical integration of I.V.P. (|3j) 

Now, using the results in Section 2 and Subsection 3.1, we can enounce the main result of this paper 
Theorem 4 Let I. V.P. ^ with g continuous. Then the problem can be numerically integrated. 
Proof 

The I.V.P. (|3| can be transformed into a set-valued I.V.P., via Filippovs regularization (Section [2]) 
which can be continuously approximated (Proposition pi) . The obtained single- valued I.V.P. may be 



numerically integrated via DFF (Section 3.1) 



One of the most unpleasant aspects of the DFF scheme is that it requires each step, the entire 
backward integration history. In other words, each value Xk is a determined function of all previous 
values xq, xi, . . . , x^-i 

or, for a large k, this implies supplementary computational cost. Thus, unlike the classical deriva- 
tives of integer order, the fractional derivatives operators are not local (i.e. we cannot determine 
D'^x{t) by using only a few values of x in a neighborhood of t). 

Remark 5 Lipschitz continuity of g is necessary only to ensure the uniqueness either of the I.V.P. 
) or ^ if the study of convergence properties and errors of some numerical method is required. 

Error analysis is a difficult task, especially because of the numerical integration of fractional equations 
(see a detailed study in [21]). Therefore, in our study we have chosen empirically h so that, for the 
common case q — 1, the obtained attractor fits as well possible the known attractor shape in the 
phase space M.^. Also, the optimal choice of the step length must ensure the maximum accuracy in the 
approximate solution at minimum computational cost. 

In order to avoid the specific phenomena for discontinuous d.s. in the discontinuity surfaces (e.g. 
sliding modes [IQ]), the initial conditions need to be chosen in Di. As can be seen from the obtained 
images, the trajectories present some "corners" typical of discontinuous d.s. 



4 Examples 

Let us first introduce the following notion 

Notation Let us denote by d = 3q the order of the d.s. modeled by I.V.P. [3] and by d* the lowest 
value for which chaos persists. 

While in the case of the continuous d.s. of integer-order chaos can appear only in nonlinear systems 
with order minimum 3, in nonlinear systems of fractional-order, this is not the case. For example, in 
|33| has been shown that a Chua circuit of order d as low as 2.7 can produce a chaotic attractor; in 
|34| it is shown that a sinusoidally nonautonomous driven Duffing system with d less than 2 can still 
behave in a chaotic manner; in |35| it is proved that the lowest order d for fractional Lorenz system 
is d* = 2.97 and in j2H] it is proved that the lowest-order chaotic system out of all the found chaotic 
systems reported in specialized literature, appears in the fractional Lii system. 

In this section we explore the chaos disappearance, by searching d* with the algorithm proposed 
before, in two examples, of known discontinuous three-dimensional d.s. In the simulations, we only 
vary the fractional derivative order q while the control parameter p is fixed so the underlying d.s. of 
integer order behaves chaotically. Simulations are performed by lowering q from q — I (the integer 
case) until chaos behavior disappears. 
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(a) 



(b) 

E=o.ai 




All computations in this paper were performed in double precision arithmetic, N was chosen large 
enough (of order 10E4) to provide accurate results, the integration step size h of IE — 3 and e = IE — 3 
of order. 

Remark 6 In all the practical examples we have found in the scientific texts ( either two-dimensional, 



typical for dry friction problems, or three-dimensional, found generally in many areas of science) there 
is at most one signum function in each equation of I. V.P. This helps to facilitate the numerical 
integration (as well known, it is time-consuming to compute the exponential functions). 



4.1 Sprott system 

The first example considered in this paper was proposed by Sprott in |32j and [38j and it is a repre- 
sentative one for the considered class of d.s. 



d^xi 
dti 

d''X2 

dti 
d^ 
dti 

whith 



X2, (14) 
X3, 

-xi - X2- px3 + sgn{xi). 



X2 

g{x) = I X3 



-Xi -X2- PX3 





(0 



















^ 1 








In this case, the chaos corresponding to p = 0.5 (see Fig. [2] a for g = 1 and Fig. [2]b for g = 0.9) 
vanishes at g = 0.89 when d* = 2.57 and when an attractive fixed point is obtained (Fig. [2]c). 



4.2 Chen system 

The next example is the fractional variant of the discontinuous Chen system 'Wf] belonging to a more 
general class of discontinuous fractional systems 

-p{x2~Xi), (15) 

sgn{xi){7 - p - X3) 0.7x2, 
sgn{x2)xi — 0.168x3. 

For g = 1 the attractor (chaotic for p — 1.18) is presented in Fig. [Sja. In this chase chaos seems to 
be more "persistent" (see Fig. [3]b,c) since it exists until q ~ 0.75 (see Fig. [sjd) when a stable limit 
cycle appears. Thus, d* = 2.25. 

Conclusion In this paper, a continuous approximation, via Filippov regularization, of a class 
of fractional-order discontinuous I. V.P. has been considered. The obtained I. V.P. of fractional-order 
was integrated by using the algorithm proposed by Diethlem, Ford and Freed in |4|. Discontinuous 
fractional-order Sprott and Chen systems have been considered. We have found that chaos still exists 
in these systems if the order is lower than three. In the cases of continuous fractional systems, the 
chaotic behavior is affected when the fractional-order decreases. However, in our opinion, choosing the 
optimal value for q (and implicitly for d), so that the mathematical model ([s]) can better illustrate the 
behavior and properties of a physical system, is much more important than finding the lowest d* (a 
simple reason could be that for d = d* chaos seems to vanish, in possible contradiction with the real 
system) . 



dixi 
dti 

d''X2 

dti 

d1X3 

dti 
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